% T18 station
load /home/mnair/projects/neesha/t18.mat
load /home/mnair/projects/neesha/Fcd18.mat
% NWP station
load /nfs/satmag_work/mnair/projects/ocean/OBEM/Toh/time.mat
load /nfs/satmag_work/mnair/projects/ocean/OBEM/Toh/Fcd.mat
time(1:3) = [];
% T16 station
load /nfs/satmag_work/mnair/projects/ocean/OBEM/Baba/Fcd16.mat
load /nfs/satmag_work/mnair/projects/ocean/OBEM/Baba/t16.mat
% T13 station
load /nfs/satmag_work/mnair/projects/ocean/OBEM/Baba/Fcd13.mat
load /nfs/satmag_work/mnair/projects/ocean/OBEM/Baba/t13.mat
%T14 station
load /nfs/satmag_work/mnair/projects/ocean/OBEM/Baba/t14.mat
load /nfs/satmag_work/mnair/projects/ocean/OBEM/Baba/Fcd14.mat

%T15 stations
load /Users/manojnair/data/obem/Fcd15.mat
load /Users/manojnair/data/obem/t15.mat

%%
% periods=[1:.1:11 11:.04:13 13:.1:14 14:.04:16 16:.1:23 23:.04:25 25:.1:30];
% periods=[1:.1:30];
periods = [4.0000    4.8000    6.0000    8.0000   11.9672   12.0000   12.4210   12.6583   23.9345   24.0000   25.8910];
%
%lon= 159.963; % The longitude of the NWP OBEM
%lon = 147.1715; % t-18
%lon = 143.96; % t16
%lon = 139.3; % t13
%lon = 139.5;%T14
lon = 141.32; %T15

t = t15 + (lon/15)/24; %t15
%t = t14 + (lon/15)/24; %t14
%t = t13 + (lon/15)/24; %OK
%t = t16 + (lon/15)/24; %OK
%t = t18 + (lon/15)/24; %OK
%t = time + (lon/15)/24; %OK NWP
tf = [1:1:length(t)]; %OK

% X Matrix

x=[cos(2*pi*t(:)/(periods(1)/24) ) sin(2*pi*t(:)/(periods(1)/24) )];

for j=2:length(periods)
    x=[x cos(2*pi*t(:)/(periods(j)/24) ) sin(2*pi*t(:)/(periods(j)/24) )];
end

% 3. Night and Day separation for Fcalc

% Get indices for night time data
N=(t-floor(t)>18/24) | (t-floor(t)<6/24);



% Get indices for day time data
D=(t-floor(t)<18/24 & (t-floor(t)>6/24))';


%% Synthetic data

synx=[cos(2*pi*t(:)/(periods(1)/24) + pi/2)];

for j=2:length(periods)
    synx=[synx cos(2*pi*t(:)/(periods(j)/24) + pi/2)];
end

synx = synx + 2* randn([length(periods),length(synx)])';
%% Testing the phase issue

amp = 1;
period = 10;
x=[cos(2*pi*t(:)/(period/24) ) sin(2*pi*t(:)/(period/24) )];

for i = 0.1:0.1:2,
    amp = i;
    Xf = amp*sin(2*pi*t/(period/24))
    + Xd;
    
    a=x\Xf;
    [b,stats] = robustfit(x,Xf);
    
    fprintf('%f %f %f %f\n',i, sqrt(a(1)^2+a(2)^2), sqrt(b(2)^2+b(3)^2), sqrt(b(2)^2+b(3)^2)- sqrt(i^2+i^2));
end;

%%
y=1;
for period=11:2:13
    for amp=.1:.1:2
        Xf = amp*sin(2*pi*tf/(period*60))+Xd';
        
        Fc_spectra_all_ls(1:end,y)=x\Xf;
        %Fc_spectra_data_all_ls(1:end,y)=abs(complex(Fc_spectra_all_ls(1:2:end,y),Fc_spectra_all_ls(2:2:end,y)));
        
        y=y+1;
    end
end

%% Plot
hold on
for y=10:20:160
    % plot(periods,Fc_spectra_data_all_ls(1:end,y))
    % plot(periods,Fc_spectra_data_day_ls(1:end,y))
    % plot(periods,Fc_spectra_data_night_ls(1:end,y))
end
axis([0 16 0 60])
xlabel('Period (hrs)')
ylabel('Amplitude (nT)')
title('Least Squares Test: Night Only')


%% Day only
Fc_spectra_day_ls=zeros(774,40);
Fc_spectra_data_day_ls=zeros(387,40);

for y=1:160
    b= Xf(1:end,y);
    Fc_spectra_day_ls(1:end,y)=x(~N,:)\b(~N);
    Fc_spectra_data_day_ls(1:end,y)=abs(complex(Fc_spectra_day_ls(1:2:end,y),Fc_spectra_day_ls(2:2:end,y)));
end

%% Night only

Fc_spectra_night_ls=zeros(6,40);
Fc_spectra_data_night_ls=zeros(3,40);

for y=1:40
    b= Xf(1:end,y);
    Fc_spectra_night_ls(1:end,y)=x(~N,:)\b(~N);
    Fc_spectra_data_night_ls(1:end,y)=abs(complex(Fc_spectra_night_ls(1:2:end,y),Fc_spectra_night_ls(2:2:end,y)));
end

%%
% Amplitudes Represented by Red *: 10, Green *: 5, Blue *:
% 4, Cyan *: 3, Magenta *: 2, Black *: 100, Red o: 1, Blue 0: 0.5

figure(5)
hold on
% plot(periods,F100_spectra_data_day_ls,'*','MarkerEdgeColor','k','MarkerSize',10)
% plot(periods,F100_spectra_data_night_ls,'*','MarkerEdgeColor','k','MarkerSize',10)

% plot(periods,F10_spectra_data_day_ls,'*','MarkerEdgeColor','r','MarkerSize',10)
% plot(periods,F10_spectra_data_night_ls,'*','MarkerEdgeColor','r','MarkerSize',10)

% plot(periods,F5_spectra_data_day_ls,'*','MarkerEdgeColor','g','MarkerSize',10)
% plot(periods,F5_spectra_data_night_ls,'*','MarkerEdgeColor','g','MarkerSize',10)

% plot(periods,F4_spectra_data_day_ls,'*','MarkerEdgeColor','b','MarkerSize',10)
% plot(periods,F4_spectra_data_night_ls,'*','MarkerEdgeColor','b','MarkerSize',10)

% plot(periods,F3_spectra_data_day_ls,'*','MarkerEdgeColor','c','MarkerSize',10)
% plot(periods,F3_spectra_data_night_ls,'*','MarkerEdgeColor','c','MarkerSize',10)

% plot(periods,F2_spectra_data_day_ls,'*','MarkerEdgeColor','m','MarkerSize',10)
% plot(periods,F2_spectra_data_night_ls,'*','MarkerEdgeColor','m','MarkerSize',10)

plot(periods,F1_spectra_data_day_ls,'*','MarkerEdgeColor','r','MarkerSize',10)
% plot(periods,F1_spectra_data_night_ls,'*','MarkerEdgeColor','r','MarkerSize',10)

plot(periods,F05_spectra_data_day_ls,'o','MarkerEdgeColor','b','MarkerSize',10)
% plot(periods,F05_spectra_data_night_ls,'o','MarkerEdgeColor','b','MarkerSize',10)

% Constant Spacing
% yday=periods*0.0856-0.0741;
% plot(periods,yday,'g','LineWidth',2)
% ynight=(0.121*periods-0.1007);
% plot(periods,ynight,'k')

% Variable Spacing
% yday=periods*0.0501-0.0439;
% plot(periods,yday,'k','LineWidth',2)
% ynight=(0.0579*periods-0.0382);
% plot(periods,ynight,'k','LineWidth',2)

axis([0 30 0 10])

xlabel('Period (hours)')
ylabel('Amplitude (nT)')
title('LSI for the test, day only.')

%% Get the scalar magnetic anomaly from the predicted vector components
% This is done by projecting the data on to the IGRF 2002 X,Y,and Z


% First get the IGRF data. This data was generated by igrf_grid.exe
% program. Each file contains one component at 0.25 x 0.25 global
% resolution. The bounds are lat -89.875 to 89.875, lon -179.875 to 179.875

igrfx = load('/nfs/satmag_work/mnair/projects/ocean/OBEM/IGRF_2002_X_Y_Z/IGRF2002X.txt');
igrfy = load('/nfs/satmag_work/mnair/projects/ocean/OBEM/IGRF_2002_X_Y_Z/IGRF2002Y.txt');
igrfz = load('/nfs/satmag_work/mnair/projects/ocean/OBEM/IGRF_2002_X_Y_Z/IGRF2002Z.txt');

Bx = reshape(igrfx(:,3),[1440,720]);
By = reshape(igrfy(:,3),[1440,720]);
Bz = reshape(igrfz(:,3),[1440,720]);

clear igrfx igrfy igrfz;

% Now rotate the matrix and cycle it to 0 to 360 interval (from -180 to
% 180)

Bz_temp = (flipud(Bz'));
Bx_temp = (flipud(Bx'));
By_temp = (flipud(By'));


% cycle rotate so that it matches with the prediction

Bz = circshift(Bz_temp,[0,720]);
Bx = circshift(Bx_temp,[0,720]);
By = circshift(By_temp,[0,720]);



% Now read the predicted data


File_1 = '/nfs/satmag_work/mnair/projects/ocean/OBEM/tpxo72_720_1440_model_results/tides_M2_720_1440.res';

fid = fopen(File_1, 'r', 'n');

LBlock = fread(fid, 1, 'long');
NO = fread(fid, 1, 'long');
Period = fread(fid, 1, 'double');
np = fread(fid, 1, 'long');
nt = fread(fid, 1, 'long');
fseek(fid, LBlock*4, 'bof');

Hpr = fread(fid, [nt np], 'float')*4*pi*100.;
Hpi = fread(fid, [nt np], 'float')*4*pi*100.;


Htr = fread(fid, [nt np], 'float')*4*pi*100.;
Hti = fread(fid, [nt np], 'float')*4*pi*100.;


Hrr = fread(fid, [nt np], 'float')*4*pi*100.;
Hri = fread(fid, [nt np], 'float')*4*pi*100.;


% now project the prediction to the main field

% create unit vector in the direction of main field

F = sqrt(Bx(:).^2+By(:).^2+Bz(:).^2);

main = [Bx(:)./F By(:)./F Bz(:)./F];

% find the scalar anomaly by projecting the predicted vector
% on to the direction of main field

% Bx = -Ht, Bz = -Hr

B_scalar_real = dot([-Htr(:) Hpr(:) -Hrr(:)]',main');
B_scalar_imag = dot([-Hti(:) Hpi(:) -Hri(:)]',main');

% reshape the matrix back

Sclar_Anly_Real =  reshape(B_scalar_real,[720,1440]);
Sclar_Anly_Imag =  reshape(B_scalar_imag,[720,1440]);


% define axes

theta = [0.5 : 1. : 180.]';
phi   = [0.5 : 1. : 360.]';

%load boundary map

load coast

% plot the map
figure(1);

imagesc(phi, theta, Sclar_Anly_Real, [-10 10] );
hold on; plot([long; 360+long], 90-[lat; lat], '-k');
set(gca, 'FontSize' , 16);
colorbar;
set(gca, 'Xtick', [0:60:360], 'YTick', [0:30:180]);
title 'Predicted scalar anomaly for M2 - real part ';
saveas(gcf,'/nfs/satmag_work/mnair/projects/ocean/OBEM/figures/Pred_M2_Scalar_Real.eps','eps');
saveas(gcf,'/nfs/satmag_work/mnair/projects/ocean/OBEM/figures/Pred_M2_Scalar_Real.fig','fig');

figure(2);

imagesc(phi, theta, Sclar_Anly_Imag, [-10 10] );
hold on; plot([long; 360+long], 90-[lat; lat], '-k');
set(gca, 'FontSize' , 16);
colorbar;
set(gca, 'Xtick', [0:60:360], 'YTick', [0:30:180]);
title 'Predicted scalar anomaly for M2 - imaginary part ';
saveas(gcf,'/nfs/satmag_work/mnair/projects/ocean/OBEM/figures/Pred_M2_Scalar_Imag.eps','eps');
saveas(gcf,'/nfs/satmag_work/mnair/projects/ocean/OBEM/figures/Pred_M2_Scalar_Imag.fig','fig');

figure(3)
%plot the scalar anomaly amplitude
imagesc(phi, theta, sqrt(Sclar_Anly_Real.^2+Sclar_Anly_Imag.^2), [0 5] )
hold on; plot([long; 360+long], 90-[lat; lat], 'k-','MarkerSize',1);
set(gca, 'FontSize' , 16);
colorbar;
set(gca, 'Xtick', [0:60:360], 'YTick', [0:30:180]);
caxis([0,5]);

saveas(gcf,'/nfs/satmag_work/mnair/projects/ocean/OBEM/figures/Pred_M2_Scalar_Amp.eps','eps');
saveas(gcf,'/nfs/satmag_work/mnair/projects/ocean/OBEM/figures/Pred_M2_Scalar_Amp.fig','fig');


save /nfs/satmag_work/mnair/projects/ocean/OBEM/tpxo72_720_1440_model_results/Scalar_M2_Anomaly_Map Sclar_Anly_Imag Sclar_Anly_Real


% Same as the above except that it is N2
%% Get the scalar magnetic anomaly from the predicted vector components
% This is done by projecting the data on to the IGRF 2002 X,Y,and Z


% First get the IGRF data. This data was generated by igrf_grid.exe
% program. Each file contains one component at 0.25 x 0.25 global
% resolution. The bounds are lat -89.875 to 89.875, lon -179.875 to 179.875

igrfx = load('/nfs/satmag_work/mnair/projects/ocean/OBEM/IGRF_2002_X_Y_Z/IGRF2002X.txt');
igrfy = load('/nfs/satmag_work/mnair/projects/ocean/OBEM/IGRF_2002_X_Y_Z/IGRF2002Y.txt');
igrfz = load('/nfs/satmag_work/mnair/projects/ocean/OBEM/IGRF_2002_X_Y_Z/IGRF2002Z.txt');

Bx = reshape(igrfx(:,3),[1440,720]);
By = reshape(igrfy(:,3),[1440,720]);
Bz = reshape(igrfz(:,3),[1440,720]);

clear igrfx igrfy igrfz;

% Now rotate the matrix and cycle it to 0 to 360 interval (from -180 to
% 180)

Bz_temp = (flipud(Bz'));
Bx_temp = (flipud(Bx'));
By_temp = (flipud(By'));


% cycle rotate so that it matches with the prediction

Bz = circshift(Bz_temp,[0,720]);
Bx = circshift(Bx_temp,[0,720]);
By = circshift(By_temp,[0,720]);



% Now read the predicted data


File_1 = '/nfs/satmag_work/mnair/projects/ocean/OBEM/tpxo72_720_1440_model_results/tides_N2_720_1440.res';

fid = fopen(File_1, 'r', 'n');

LBlock = fread(fid, 1, 'long');
NO = fread(fid, 1, 'long');
Period = fread(fid, 1, 'double');
np = fread(fid, 1, 'long');
nt = fread(fid, 1, 'long');
fseek(fid, LBlock*4, 'bof');

Hpr = fread(fid, [nt np], 'float')*4*pi*100.;
Hpi = fread(fid, [nt np], 'float')*4*pi*100.;


Htr = fread(fid, [nt np], 'float')*4*pi*100.;
Hti = fread(fid, [nt np], 'float')*4*pi*100.;


Hrr = fread(fid, [nt np], 'float')*4*pi*100.;
Hri = fread(fid, [nt np], 'float')*4*pi*100.;


% now project the prediction to the main field

% create unit vector in the direction of main field

F = sqrt(Bx(:).^2+By(:).^2+Bz(:).^2);

main = [Bx(:)./F By(:)./F Bz(:)./F];

% find the scalar anomaly by projecting the predicted vector
% on to the direction of main field

% Bx = -Ht, Bz = -Hr

B_scalar_real = dot([-Htr(:) Hpr(:) -Hrr(:)]',main');
B_scalar_imag = dot([-Hti(:) Hpi(:) -Hri(:)]',main');

% reshape the matrix back

Sclar_Anly_Real =  reshape(B_scalar_real,[720,1440]);
Sclar_Anly_Imag =  reshape(B_scalar_imag,[720,1440]);


% define axes

theta = [0.5 : 1. : 180.]';
phi   = [0.5 : 1. : 360.]';

%load boundary map

load coast

% plot the map
figure(1);

imagesc(phi, theta, Sclar_Anly_Real, [-1 1] );
hold on; plot([long; 360+long], 90-[lat; lat], '-k');
set(gca, 'FontSize' , 16);
colorbar;
set(gca, 'Xtick', [0:60:360], 'YTick', [0:30:180]);
title 'Predicted scalar anomaly for N2 - real part ';
saveas(gcf,'/nfs/satmag_work/mnair/projects/ocean/OBEM/figures/Pred_N2_Scalar_Real.eps','eps');
saveas(gcf,'/nfs/satmag_work/mnair/projects/ocean/OBEM/figures/Pred_N2_Scalar_Real.fig','fig');

figure(2);

imagesc(phi, theta, Sclar_Anly_Imag, [-1 1] );
hold on; plot([long; 360+long], 90-[lat; lat], '-k');
set(gca, 'FontSize' , 16);
colorbar;
set(gca, 'Xtick', [0:60:360], 'YTick', [0:30:180]);
title 'Predicted scalar anomaly for N2 - imaginary part ';
saveas(gcf,'/nfs/satmag_work/mnair/projects/ocean/OBEM/figures/Pred_N2_Scalar_Imag.eps','eps');
saveas(gcf,'/nfs/satmag_work/mnair/projects/ocean/OBEM/figures/Pred_N2_Scalar_Imag.fig','fig');

figure(3)
%plot the scalar anomaly amplitude
imagesc(phi, theta, sqrt(Sclar_Anly_Real.^2+Sclar_Anly_Imag.^2), [0 1] )
hold on; plot([long; 360+long], 90-[lat; lat], 'k-','MarkerSize',1);
set(gca, 'FontSize' , 16);
colorbar;
set(gca, 'Xtick', [0:60:360], 'YTick', [0:30:180]);
caxis([0,1]);
title 'Predicted scalar Amplitude for N2 '
saveas(gcf,'/nfs/satmag_work/mnair/projects/ocean/OBEM/figures/Pred_N2_Scalar_Amp.eps','eps');
saveas(gcf,'/nfs/satmag_work/mnair/projects/ocean/OBEM/figures/Pred_N2_Scalar_Amp.fig','fig');


save /nfs/satmag_work/mnair/projects/ocean/OBEM/tpxo72_720_1440_model_results/Scalar_N2_Anomaly_Map Sclar_Anly_Imag Sclar_Anly_Real


% Same as the above except that it is O1
%% Get the scalar magnetic anomaly from the predicted vector components
% This is done by projecting the data on to the IGRF 2002 X,Y,and Z


% First get the IGRF data. This data was generated by igrf_grid.exe
% program. Each file contains one component at 0.25 x 0.25 global
% resolution. The bounds are lat -89.875 to 89.875, lon -179.875 to 179.875

igrfx = load('/nfs/satmag_work/mnair/projects/ocean/OBEM/IGRF_2002_X_Y_Z/IGRF2002X.txt');
igrfy = load('/nfs/satmag_work/mnair/projects/ocean/OBEM/IGRF_2002_X_Y_Z/IGRF2002Y.txt');
igrfz = load('/nfs/satmag_work/mnair/projects/ocean/OBEM/IGRF_2002_X_Y_Z/IGRF2002Z.txt');

Bx = reshape(igrfx(:,3),[1440,720]);
By = reshape(igrfy(:,3),[1440,720]);
Bz = reshape(igrfz(:,3),[1440,720]);

clear igrfx igrfy igrfz;

% Now rotate the matrix and cycle it to 0 to 360 interval (from -180 to
% 180)

Bz_temp = (flipud(Bz'));
Bx_temp = (flipud(Bx'));
By_temp = (flipud(By'));


% cycle rotate so that it matches with the prediction

Bz = circshift(Bz_temp,[0,720]);
Bx = circshift(Bx_temp,[0,720]);
By = circshift(By_temp,[0,720]);



%% Now read the predicted data


%File_1 = '/nfs/satmag_work/mnair/projects/ocean/OBEM/tpxo72_720_1440_model_results/tides_N2_720_1440.res';
%File_1 = '/nfs/satmag_work/mnair/projects/ocean/OBEM/tpxo72_720_1440_modeling/Global_1D_Kuvshinov_Olsen_2006/Kuvsh_tides_N2_720_1440.res';
%File_1 = '/nfs/satmag_work/mnair/projects/ocean/OBEM/tpxo72_720_1440_modeling/North_Pacific_1D_section_Baba_et_al_2010/Baba_NP_tides_O1_720_1440.res';
File_1 = '/nfs/satmag_work/mnair/projects/ocean/OBEM/tpxo72_720_1440_modeling/Philippine_Sea_1D_section_Baba_et_al_2010/Baba_PS_tides_M2_720_1440.res';

fid = fopen(File_1, 'r', 'n');

LBlock = fread(fid, 1, 'long');
NO = fread(fid, 1, 'long');
Period = fread(fid, 1, 'double');
np = fread(fid, 1, 'long');
nt = fread(fid, 1, 'long');
fseek(fid, LBlock*4, 'bof');

Hpr = fread(fid, [nt np], 'float')*4*pi*100.;
Hpi = fread(fid, [nt np], 'float')*4*pi*100.;


Htr = fread(fid, [nt np], 'float')*4*pi*100.;
Hti = fread(fid, [nt np], 'float')*4*pi*100.;


Hrr = fread(fid, [nt np], 'float')*4*pi*100.;
Hri = fread(fid, [nt np], 'float')*4*pi*100.;


% now project the prediction to the main field

% create unit vector in the direction of main field

F = sqrt(Bx(:).^2+By(:).^2+Bz(:).^2);

main = [Bx(:)./F By(:)./F Bz(:)./F];

% find the scalar anomaly by projecting the predicted vector
% on to the direction of main field

% Bx = -Ht, Bz = -Hr

B_scalar_real = dot([-Htr(:) Hpr(:) -Hrr(:)]',main');
B_scalar_imag = dot([-Hti(:) Hpi(:) -Hri(:)]',main');

% reshape the matrix back

Sclar_Anly_Real =  reshape(B_scalar_real,[720,1440]);
Sclar_Anly_Imag =  reshape(B_scalar_imag,[720,1440]);

%% Plot
% define axes

theta = [0.5 : 1. : 180.]';
phi   = [0.5 : 1. : 360.]';

%load boundary map

load coast

% plot the map
figure(1);

imagesc(phi, theta, Sclar_Anly_Real, [-1 1] );
hold on; plot([long; 360+long], 90-[lat; lat], '-k');
set(gca, 'FontSize' , 16);
colorbar;
set(gca, 'Xtick', [0:60:360], 'YTick', [0:30:180]);
title 'Predicted scalar anomaly for O1 - real part ';
saveas(gcf,'/nfs/satmag_work/mnair/projects/ocean/OBEM/figures/Baba_PS_1D_720_1440_O1_Scalar_Real.eps','eps');
saveas(gcf,'/nfs/satmag_work/mnair/projects/ocean/OBEM/figures/Baba_PS_1D_720_1440_O1_Scalar_Real.fig','fig');

figure(2);

imagesc(phi, theta, Sclar_Anly_Imag, [-1 1] );
hold on; plot([long; 360+long], 90-[lat; lat], '-k');
set(gca, 'FontSize' , 16);
colorbar;
set(gca, 'Xtick', [0:60:360], 'YTick', [0:30:180]);
title 'Predicted scalar anomaly for O1 - imaginary part ';
saveas(gcf,'/nfs/satmag_work/mnair/projects/ocean/OBEM/figures/Baba_PS_1D_720_1440_O1_Scalar_Imag.eps','eps');
saveas(gcf,'/nfs/satmag_work/mnair/projects/ocean/OBEM/figures/Baba_PS_1D_720_1440_O1_Scalar_Imag.fig','fig');

figure(3)
%plot the scalar anomaly amplitude
imagesc(phi, theta, sqrt(Sclar_Anly_Real.^2+Sclar_Anly_Imag.^2), [0 5] )
hold on; plot([long; 360+long], 90-[lat; lat], 'k-','MarkerSize',1);
set(gca, 'FontSize' , 16);
colorbar;
set(gca, 'Xtick', [0:60:360], 'YTick', [0:30:180]);
caxis([0,1]);
title 'Predicted scalar Amplitude for O1 '
saveas(gcf,'/nfs/satmag_work/mnair/projects/ocean/OBEM/figures/Baba_PS_1D_720_1440_O1_Scalar_Amp.eps','eps');
saveas(gcf,'/nfs/satmag_work/mnair/projects/ocean/OBEM/figures/Baba_PS_1D_720_1440_O1_Scalar_Amp.fig','fig');


save /nfs/satmag_work/mnair/projects/ocean/OBEM/tpxo72_720_1440_model_results/Baba_PS_1D_720_1440_O1_Scalar_Anomaly_Map Sclar_Anly_Imag Sclar_Anly_Real


%% Adding contour of Bz as per reviewer's suggestion

uiopen('/nfs/satmag_work/mnair/projects/ocean/OBEM/figures/Pred_M2_Scalar_Amp.fig',1);
[X_1440,Y_1440] = meshgrid(0.125:0.25:360,0.125:0.25:180);
[cs, h ] = contour(X_1440, Y_1440, Bz,'color','black');
clabel(cs, h,'FontSize',16);

%% Get the scalar magnetic anomaly from the predicted vector components180x360 modeling
% This is done by projecting the data on to the IGRF 2002 X,Y,and Z


% First get the IGRF data. This data was generated by igrf_grid.exe
% program. Each file contains one component at 0.25 x 0.25 global
% resolution. The bounds are lat -89.875 to 89.875, lon -179.875 to 179.875

igrfx = load('/nfs/satmag_work/mnair/projects/ocean/OBEM/IGRF_2002_X_Y_Z/IGRF2002X.txt');
igrfy = load('/nfs/satmag_work/mnair/projects/ocean/OBEM/IGRF_2002_X_Y_Z/IGRF2002Y.txt');
igrfz = load('/nfs/satmag_work/mnair/projects/ocean/OBEM/IGRF_2002_X_Y_Z/IGRF2002Z.txt');

Bx = reshape(igrfx(:,3),[1440,720]);
By = reshape(igrfy(:,3),[1440,720]);
Bz = reshape(igrfz(:,3),[1440,720]);

clear igrfx igrfy igrfz;

% Now rotate the matrix and cycle it to 0 to 360 interval (from -180 to
% 180)

Bz_temp = (flipud(Bz'));
Bx_temp = (flipud(Bx'));
By_temp = (flipud(By'));


% cycle rotate so that it matches with the prediction

Bz = circshift(Bz_temp,[0,720]);
Bx = circshift(Bx_temp,[0,720]);
By = circshift(By_temp,[0,720]);


[X_360,Y_360] = meshgrid(0.5:1:360,(-90+0.5):1:90);

[X_1440,Y_1440] = meshgrid(0.125:0.25:360,(-90+0.125):0.25:90);


% get the grid at 180x360

Bz_360 = interp2(X_1440,Y_1440,Bz,X_360,Y_360);
By_360 = interp2(X_1440,Y_1440,By,X_360,Y_360);
Bx_360 = interp2(X_1440,Y_1440,Bx,X_360,Y_360);



%% Now read the predicted data
%clear;
% Sclar_Anly_Imag_new = Sclar_Anly_Imag;
% Sclar_Anly_Real_new = Sclar_Anly_Real;

%File_1 = '/nfs/satmag_work/mnair/projects/ocean/OBEM/tpxo72_180_360_modeling/sh_1D_Shimizu_lateral_lateral_lith/tides_M2_180_360_3K_300_slab.res';
%File_1 = '/nfs/satmag_work/mnair/projects/ocean/OBEM/tpxo72_180_360_modeling/sh_1D_Shimizu_lateral_lateral_lith/tides_M2_180_360_3K_300.res';
%File_1 = '/nfs/satmag_work/mnair/projects/ocean/OBEM/tpxo72_180_360_modeling/sh_1D_Shimizu_lateral_lateral_lith/tides_O1_180_360.res';
%File_1 = '/nfs/satmag_work/mnair/projects/ocean/OBEM/tpxo72_180_360_modeling/sh_1D_Shimizu_lateral_lateral_lith/tides_O1_180_360_30K_3K.res';
File_1 = '/nfs/satmag_work/mnair/projects/ocean/OBEM/tpxo72_180_360_modeling/sh_1D_Shimizu/tides_O1_180_360.res';


fid = fopen(File_1, 'r', 'n');

LBlock = fread(fid, 1, 'long');
NO = fread(fid, 1, 'long');
Period = fread(fid, 1, 'double');
np = fread(fid, 1, 'long');
nt = fread(fid, 1, 'long');
fseek(fid, LBlock*4, 'bof');

Hpr = fread(fid, [nt np], 'float')*4*pi*100.;
Hpi = fread(fid, [nt np], 'float')*4*pi*100.;


Htr = fread(fid, [nt np], 'float')*4*pi*100.;
Hti = fread(fid, [nt np], 'float')*4*pi*100.;


Hrr = fread(fid, [nt np], 'float')*4*pi*100.;
Hri = fread(fid, [nt np], 'float')*4*pi*100.;


% now project the prediction to the main field

% create unit vector in the direction of main field

F = sqrt(Bx_360(:).^2+By_360(:).^2+Bz_360(:).^2);

main = [Bx_360(:)./F By_360(:)./F Bz_360(:)./F];

% find the scalar anomaly by projecting the predicted vector
% on to the direction of main field

% Bx = -Ht, Bz = -Hr

B_scalar_real = dot([-Htr(:) Hpr(:) -Hrr(:)]',main');
B_scalar_imag = dot([-Hti(:) Hpi(:) -Hri(:)]',main');

% reshape the matrix back

Sclar_Anly_Real =  reshape(B_scalar_real,[180,360]);
Sclar_Anly_Imag =  reshape(B_scalar_imag,[180,360]);


%% define axes

theta = [0.5 : 1. : 180.]';
phi   = [0.5 : 1. : 360.]';

%load boundary map

load coast

% plot the map
figure(1);

imagesc(phi, theta, Sclar_Anly_Real, [-5 1] );
hold on; plot([long; 360+long], 90-[lat; lat], '-k');
set(gca, 'FontSize' , 16);
colorbar;
set(gca, 'Xtick', [0:60:360], 'YTick', [0:30:180]);
title 'Predicted scalar anomaly for M2(180x360) - real part ';
saveas(gcf,'/nfs/satmag_work/mnair/projects/ocean/OBEM/figures/Pred_180_360_M2_Scalar_Real.eps','eps');
saveas(gcf,'/nfs/satmag_work/mnair/projects/ocean/OBEM/figures/Pred_180_360_M2_Scalar_Real.fig','fig');

figure(2);

imagesc(phi, theta, Sclar_Anly_Imag, [-1 1] );
hold on; plot([long; 360+long], 90-[lat; lat], '-k');
set(gca, 'FontSize' , 16);
colorbar;
set(gca, 'Xtick', [0:60:360], 'YTick', [0:30:180]);
title 'Predicted scalar anomaly for M2(180x360) - imaginary part ';
saveas(gcf,'/nfs/satmag_work/mnair/projects/ocean/OBEM/figures/Pred_180_360_M2_Scalar_Imag.eps','eps');
saveas(gcf,'/nfs/satmag_work/mnair/projects/ocean/OBEM/figures/Pred_180_360_M2_Scalar_Imag.fig','fig');

figure(3)
%plot the scalar anomaly amplitude
imagesc(phi, theta, sqrt(Sclar_Anly_Real.^2+Sclar_Anly_Imag.^2), [0 5] )
hold on; plot([long; 360+long], 90-[lat; lat], 'k-','MarkerSize',1);
set(gca, 'FontSize' , 16);
colorbar;
set(gca, 'Xtick', [0:60:360], 'YTick', [0:30:180]);
caxis([0,5]);
title 'Predicted scalar Amplitude for M2(180x360) '
saveas(gcf,'/nfs/satmag_work/mnair/projects/ocean/OBEM/figures/Pred_180_360_M2_Scalar_Amp.eps','eps');
saveas(gcf,'/nfs/satmag_work/mnair/projects/ocean/OBEM/figures/Pred_180_360_M2_Scalar_Amp.fig','fig');


save /nfs/satmag_work/mnair/projects/ocean/OBEM/tpxo72_180_360_modeling_results/Scalar_180_360_M2_Anomaly_Map Sclar_Anly_Imag Sclar_Anly_Real






%% Station locations
station = {'NWP' 41.10 159.96 5580 08/01/01  08/01/02 90 ;...
    'T13' 24.98 139.30 4794 10/01/05  11/30/06 90 ;...
    'T14' 22.00 139.50 4945 10/01/05  11/30/06 90 ;...
    'T15' 29.00 141.32 4026 10/01/05  11/30/06 90 ;...
    'T16' 32.52 143.96 5408 11/01/07  11/30/08 90 ;...
    'T18' 27.14 147.17 5594 11/01/06  11/30/07 90 };


station_lat = [ 41.10  24.98 22.00 29.00 32.52 27.14 ];
station_long = [ 159.96 139.30 139.50  141.32 143.96  147.17 ];
station_names = ['NWP';'T13';'T14';'T15';'T16';'T18'];



%% Create figure for tidal flow on 0.25x0.25 and 1x1 degree resolution

% Get M2 data


ncid = netcdf.open('/nfs/satmag_work/mnair/projects/ocean/OBEM/tpxo72_model_data/DATA/u_tpxo7.2.nc','NOWRITE');

% Variables
% 1 lon_u
% 2 lat_u
% 3 lon_v
% 4 lat_v
% 5 Ua
% 6 ua
% 7 up
% 8 Va
% 9 va
% 10 vp
% 11 URe
% 12 UIm
% 13 VRe
% 14 VIm

lon_u = netcdf.getVar(ncid, 1);
lat_u = netcdf.getVar(ncid, 2);
lon_v = netcdf.getVar(ncid, 3);
lat_v = netcdf.getVar(ncid, 4);
URe = netcdf.getVar(ncid, 11);
UIm = netcdf.getVar(ncid, 12);
VRe = netcdf.getVar(ncid, 13);
VIm = netcdf.getVar(ncid, 14);

% interploate the data to a common grid (the original grid is staggered)

latmat = lat_v;
lonmat = lon_u;

% remove the -90.125 values from the latmat

latmat(1,:) = [];
lonmat(1,:) = [];

% Dimention 1 is M2
% 1 M2,2 S2,3 N2,4 K2,5 K1,6 O1,7 P1,8 Q1
tidal_comp = ['M2';'S2';'N2';'K2';'K1';'O1';'P1';'Q1'];
kk = 1;


URe_inter = interp2(lon_u,lat_u, squeeze(URe(:,:,kk)), lonmat,latmat);
UIm_inter = interp2(lon_u,lat_u, squeeze(UIm(:,:,kk)), lonmat,latmat);
VRe_inter = interp2(lon_v,lat_v, squeeze(VRe(:,:,kk)), lonmat,latmat);
VIm_inter = interp2(lon_v,lat_v, squeeze(VIm(:,:,kk)), lonmat,latmat);

% select region of interest

L = lonmat > 130 & lonmat < 170 & latmat > 18 & latmat < 46;

lonnew = reshape(lonmat(L),[112,160]);
latnew = reshape(latmat(L),[112,160]);
iurenew = reshape(URe_inter(L),[112,160]);
iuimnew = reshape(UIm_inter(L),[112,160]);
ivrenew = reshape(VRe_inter(L),[112,160]);
ivimnew = reshape(VIm_inter(L),[112,160]);


[X_360,Y_360] = meshgrid(0.5:1:360,(-90+0.5):1:90);
URe_inter_1x1 = interp2(lonmat,latmat,URe_inter,X_360,Y_360);
UIm_inter_1x1 = interp2(lonmat,latmat,UIm_inter,X_360,Y_360);
VRe_inter_1x1 = interp2(lonmat,latmat,VRe_inter,X_360,Y_360);
VIm_inter_1x1 = interp2(lonmat,latmat,VIm_inter,X_360,Y_360);

% select region of interest

L = X_360 > 130 & X_360 < 170 & Y_360 > 18 & Y_360 < 46;

lonnew_1x1 = reshape(X_360(L),[28,40]);
latnew_1x1 = reshape(Y_360(L),[28,40]);
iurenew_1x1 = reshape(URe_inter_1x1(L),[28,40]);
iuimnew_1x1 = reshape(UIm_inter_1x1(L),[28,40]);
ivrenew_1x1 = reshape(VRe_inter_1x1(L),[28,40]);
ivimnew_1x1 = reshape(VIm_inter_1x1(L),[28,40]);

%% plot the Amplitude of U and V in 0.25x0.25 degree resolution


load coast

subplot(221);
cla;set(gca, 'FontSize',16);
[cs, h ] = contour(lonnew, latnew, sqrt(iurenew.^2 + iuimnew.^2),  [10:20:140],'LineWidth',2);
hold on;
plot(station_long, station_lat,'k.','MarkerSize',10);
text(station_long, station_lat,station_names,'FontSize',16);
caxis([0,150]);
h=colorbar('FontSize',16);
ylabel(h, 'm^2/s')
title('Amplitude of M2 flow (U) in 0.25x0.25 deg grid','FontSize',16);
plot(long,lat,'k','LineWidth',2);
subplot(223);
cla;
set(gca, 'FontSize',16);
[cs, h ] = contour(lonnew, latnew, sqrt(ivrenew.^2 + ivimnew.^2),  [10:20:140],'LineWidth',2);
hold on;  caxis([0,150]);
plot(station_long, station_lat,'k.','MarkerSize',10);
text(station_long, station_lat,station_names,'FontSize',16);
h=colorbar('FontSize',16);
ylabel(h, 'm^2/s')
title('Amplitude of M2 flow (V) in 0.25x0.25 deg grid','FontSize',16);
plot(long,lat,'k','LineWidth',2);

% make 1x1 grid

% plot the Amplitude of U and V in 1x1 degree resolution

subplot(222);

cla;
set(gca, 'FontSize',16);
[cs, h ] = contour(lonnew_1x1, latnew_1x1, sqrt(iurenew_1x1.^2 + iuimnew_1x1.^2),  [10:20:140],'LineWidth',2);
hold on;  caxis([0,150]);
plot(station_long, station_lat,'k.','MarkerSize',10);
text(station_long, station_lat,station_names,'FontSize',16);
h=colorbar('FontSize',16);
ylabel(h, 'm^2/s')
title('Amplitude of M2 flow (U) in 1x1 deg grid','FontSize',16);
plot(long,lat,'k','LineWidth',2);

subplot(224);
cla;
set(gca, 'FontSize',16);
[cs, h ] = contour(lonnew_1x1, latnew_1x1, sqrt(ivrenew_1x1.^2 + ivimnew_1x1.^2),  [10:20:140],'LineWidth',2);
hold on;  caxis([0,150]);
plot(station_long, station_lat,'k.','MarkerSize',10);
title('Amplitude of M2 flow (V) in 1x1 deg grid','FontSize',16);

text(station_long, station_lat,station_names,'FontSize',16);
h=colorbar('FontSize',16);
ylabel(h, 'm^2/s')
plot(long,lat,'k','LineWidth',2);

